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Abstract 



We present results of numerical studies of the two dimensional XY model 
with four and eight fold symmetry breaking fields. This model has recently 
been shown to describe hydrogen induced reconstruction on the VF(IOO) sur- 
face. Based on mean-field and renormalization group arguments, we first 
show how the interplay between the anisotropy fields can give rise to differ- 
ent phase transitions in the model. When the fields are compatible with each 
other there is a continuous phase transition when the fourth order field is var- 
ied from negative to positive values. This transition becomes discontinuous 
at low temperatures. These two regimes are separated by a multicritical 
point. In the case of competing four and eight fold fields, the first order 
transition at low temperatures opens up into two Ising transitions. We then 
use numerical methods to accurately locate the position of the multicritical 
point, and to verify the nature of the transitions. The different techniques 
used include Monte Carlo histogram methods combined with finite size scal- 
ing analysis, the real space Monte Carlo Renormalization Group method, 
and the Monte Carlo Transfer Matrix method. Our numerical results are in 
good agreement with the theoretical arguments. 



PACS numbers: 64.60. Cn, 68.35.Rh, 02.70.Lq. 



1 Introduction 



In two dimensions, conventional long range order cannot exist in continuous 
spin models (0(n),n > 2) because it is destroyed by spin wave excitations 
P], H ]3[]. However, Kosterlitz and Thouless [[!]] proposed that in the XY 
(0(2)) model there is a phase below the critical temperature where topolog- 
ical long range order can be defined. The vanishing of this order occurs via 
the Kosterlitz-Thouless (KT) transition. Physical systems where the KT 
transition occurs are numerous; they include superfluid 4 He films, Joseph- 
son junction arrays, superconducting transition of type II, and various phase 
transitions on surfaces and adsorption layers. 

In many cases, the realization of the XY model is accompanied by various 
symmetry breaking fields, whose effect is very complicated as demonstrated 
qualitatively by Jose et al. || . For example, it was recognized already at an 
early stage || that the presence of a four fold field restores a conventional 
phase transition, but with continuously varying critical exponents. In con- 
trast, a six fold symmetry breaking field opens up the KT transition into 
two parts: at high temperatures, the transition remains XY type, but at low 
enough temperatures it is into a discrete planar phase in which the system 
orders along one of the six preferred directions ||. Both of these situations 
have been realised in experimental systems; four fold fields are know to be 
present in surface structural phase transitions whereas liquid crystals provide 
a case where the six fold field exists. 
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The effects of symmetry breaking fields are further complicated by the exis- 
tence of higher order multiples of the fields, which are allowed by symmetry. 
In most cases, since these fields are irrelevant if lower order fields are present, 
their influence has been neglected. However, if it happens that the lowest 
order symmetry breaking field vanishes, the higher harmonics can become 
relevant at low enough temperatures M. A demonstration of this fact is the 
study of Selinger and Nelson || who modeled a phase transition occuring in 
liquid crystals by an XY model with six and twelve fold symmetry break- 
ing fields. They found a rich behavior of the phase diagram depending on 
whether or not the two fields are compatible with each other. 

We have recently developed |J a lattice Hamiltonian for the adsorption sys- 
tem H/W (100), which is based on the XY model with a four fold symmetry 
breaking field, and its higher harmonics. We have argued that the essential 
physics of this model is dictated by an interplay between the four fold and 
the eight fold fields, in a manner very similar to that of Ref. p|. The in- 
triguing aspect of this model is that the strengths of the symmetry fields - 
the four fold field in particular — are tunable by changing the amount of 
adsorbed hydrogen. In fact, it was demonstrated that the fourth order field 
vanishes at a hydrogen coverage of about 0.1 for this system. This system 
thus provides an ideal example for studying the effect of interplay between 
symmetry breaking fields within the XY model. 

The purpose of the present work is to conduct a detailed, quantitative study 
of the two dimensional XY model with four and eight fold anisotropy fields. 
We shall first discuss in detail, how the interplay between these two anisotropy 
fields dictates the nature of the phase diagram at low temperatures where 
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the eight fold field is relevant in the renormalization group sense. We give 
both mean field and renormalization group arguments in explaining how it 
is possible to obtain either a discontinuous, or two continuous Ising transi- 
tions at low temperatures due to the interplay between the anisotropy fields. 
Namely, when the four and eight fold fields are compatible with each other 
and do not compete, this gives rise to a first order transition at h± = as we 
pass from negative (positive) values of the four fold field to positive (nega- 
tive) with a finite hs field. However, when the two fields are competing with 
each other, the first order transition opens up into two Ising transitions with 
an intermediate phase in between. 

Following analytic arguments, we proceed to simulate the XY model with 
symmetry breaking fields using the Monte Carlo method with the Wolff 
updating algorithm which is generalized to include contributions from the 
anisotropy fields. We employ finite size scaling arguments to locate the mul- 
ticritical point in this model. These results are further corroborated by Monte 
Carlo Transfer Matrix (MCTM) studies. We also verify the existence of the 
continuous Ising transition in the case of competing anisotropy fields by using 
both the MCTM and real space Monte Carlo Renormalization Group meth- 
ods. Our results are in good agreement with theoretical predictions, and also 
consistent with available experimental data for the H/W (100) adsorption 
system. 
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2 Model Hamiltonian and Qualitative Renor- 
malization Group Analysis 



The Hamiltonian of the XY model with four and eight fold anisotropy fields 
can be written as 

H = -K c os(0i - 4>j) - h^Y^ cos(40j) + h 8 ^ cos(80;), (1) 

<i,j> i i 

where K is the XY coupling constant between nearest neighbors, <pi are the 
angle variables defined by the individual spin vectors = (cos 0,, sin 0j) on 
site i, hi and h s are the four and eight fold symmetry breaking fields, respec- 
tively, and we have subsumed the temperature into the coupling constants 
and fields. The summation < i,j > goes over the nearest neighbors, and the 
summations i are over all lattice sites. 

We will first discuss mean field theory to obtain a qualitative picture of the 
critical phenomena that Eq. (|l|) gives rise to. Our purpose is to give insight 
to the underlying physics which is dictated by the interplay between the 
anisotropy fields. These arguments correspond to the case where both fields 
are always assumed to be relevant, and will thus give correct qualitative be- 
havior at low temperatures only. There are two possible scenarios depending 
on the anisotropy potential V(<f>) of Eq. (|I|), defined by 

V(<f>) = -hi cos(40) + h 8 cos(80). (2) 

Namely, when h± is negative it favors spins aligning along the directions 
cf) = 7r/4 + nTc/2,n = 0,1,... whereas a positive h^ favors <fi = mr/2,n = 
0,1,... . The eight fold field has two possibilities as well. When h 8 < 
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the favored orientations are = rar/4, n = 0,1,... which are exactly the 
directions favored by either a negative or a positive four fold field. We say 
that in this case the anisotropy fields are non-competing or compatible with 
each other. For h$ > the favored orientations are = n / 8 + nir / 4, n = 
0,1,.... As these are different from those favored by the four fold field, 
competition will set in when the four and eight fold fields are of the same 
order of magnitude. 

We first consider the case of a negative, i.e. non-competing h s . In Fig. 1(a), 
we show the anisotropy potential V(<j)) as a function of the angle for various 
values of h 4 , given a fixed h%. The local minima at = 0, and at = ±7r/4 
are clearly visible. When h 4 > 0, the minima at = are deeper, whereas 
for h 4 < the minima at = ir/4 are deeper. As h 4 passes through zero 
from h 4 > to h 4 < 0, there is a first order transition from one minimum to 
another. 

In the case of a positive or competing eight fold field, the situation is more 
complicated. The behavior of the potential V{4>) for fixed h 8 > and various 
values of h 4 is shown in Fig. 1(b). For h 4 > 4|/i 8 |, there is only one minimum 
at = 0, and one at = ±7r/4 for h 4 < — 4|/i 8 |. But now as the four fold field 
passes through h 4 = 4|/t 8 |, the single minimum splits into two at = ±0q. 
By minimizing the anisotropy potential V(<f)), we can show that these new 
minima begin to form at 




= ±0' o . 



(3) 
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At the other boundary where hi passes through — 4|/i 8 |, the same argument 
applies except that (f)' now measures the deviation from 7r/4. The first order 
transition for the compatible field case (cf. Fig. 2(a)) has now opened up 
into two continuous transitions with an intermediate phase in between (cf. 
Fig. 2(b)). 

We will next show explicitly that the continuous phase transition belongs to 
the universality class of the two dimensional Ising model. We can approach 
the transition boundaries either from the phase where h 4 > where the 
preferred orientation of spins is along = 0, n/2, ... , or from the phase where 
hi < and the favored directions are then = ir/4, 3ir/2, ... . Let us consider 
the former possibility. For large h 4 , the system is in the energy minimum 
at = 0. When we approach the transition (h 4 — > 0), two minima form at 
= ±0q. Let us consider the contribution to the partition function due to 
these new minima. The Boltzmann weights can be written as exp{i^ cos[(sj — 
Sj)0o]} where Sj = ±1. We denote these weights by W[s] where s = Sj — Sj = 
0, 2. For two neighboring spins at the same minimum, 



W[0] = exp(K), (4) 

while for two opposite spins, 



W[2] = exp(JTcos20( ) ). (5) 

The equivalent Ising coupling constant is thus J = K sin 2 0q. Since the 
argument is valid in the limit K — > 00 only, we will later numerically verify 
the nature of the Ising transition predicted here. 
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Next, we present more quantitative renormalization group (RG) arguments 
which are a direct extension of the work by Selinger and Nelson |J. At high 
temperatures (where K is small), only the four fold field is relevant in the 
RG sense, and the long range order in the system is dictated by h± alone. At 
some temperature corresponding to K m the eight fold field becomes relevant. 
For lower temperatures (larger values of K) and for finite and hg, the 
nature of the phase diagram is determined by the interplay between these 
two fields as qualitatively discussed above. Any anisotropy field h p of order 
p will obey the RG recursion relation [Q: 

h' p = b x <>h p , (6) 
where h' p is the field obtained after an RG iteration, b is a constant, and 

*> = *- ^ (7 > 

When the parameter X p > 0, the field h' will increase as iterations proceed, 
and the p th anisotropy field is a relevant variable However, as the tem- 
perature increases, higher order fields become irrelevant (X p < 0), and the 
respective ti will decrease with iteration. On the other hand, it can easily 
be seen that at K = oo all symmetry breaking perturbations are relevant. If 
\ p = 2 — p 2 / (4nK) = 0, we get the temperature K^ 1 below which the field 
h p is relevant, i.e. it is marginally relevant at that point [|J: 



K m]p = ^f- ( 8 ) 

We can conclude that the four and eight fold anisotropy fields are relevant 
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at all temperatures K^ 1 < = it/2 and K" 1 < K~ 8 = vr/8, respectively. 
The results for the XY model with only the four fold field are well known 
M. The continuous phase transition into an ordered state (h^ 7^ 0) along the 
critical lines terminating at K^} T < K m \ = n/2 belongs to the universality 
class of the XY model with a cubic anisotropy field, where K~^ T is the (true) 
order-disorder transition temperature for the pure XY model. Along the 
critical line /j 4 = we have a series of continuous XY transitions for all 
K^ 1 < K^ T (cf. Fig. 2). In the other limit where h± — ► oo, the model 
becomes a four state clock model which can be shown to decouple into two 
Ising models 0. 

With the inclusion of the eight fold field, different scenarios occur depending 
on the signs of and h 8 . For the remainder of this paper, we will denote 
by i^" 1 the true value of the coupling constant, where the eight fold field 
becomes relevant. For K~ l < if" 1 , both and h$ are relevant. If h$ is 
positive, the renormalized potential V*(4>) is similar to that of Fig. 1(b) 
with two continuous Ising transitions. In the limit K — > 00, we obtain the 
exact mean field result: the transition occurs at = ±4|/tg|. In the other 
limit K — > K m , it can be shown || that to lowest order in \h 8 \, we have 

K^ = ^ + B\hs\, (9) 

where B > is a constant which can be estimated to be 0(1). The upper 
bound for K^ 1 is the Kosterlitz-Thouless transition temperature Kg T < n/2 
for this model. The lower bound, on the other hand, is given by the zero field 
estimate K~ 1 (/i 8 = 0) = ir/8. The corresponding phase diagram is shown 
schematically in Fig. 2(b). 
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For a negative hs and for K~ x < K^ 1 , the system fluctuates in the minima of 
Fig. 1(a). The location of the deepest minimum changes as h± passes through 
/14 = 0, and we expect a first order phase transition. Selinger and Nelson 
H] have in fact shown that in this case there is a discontinuity in the order 
parameter across the transition. This discontinuity vanishes exponentially 
as K — > K m H thus making it very difficult to numerically locate K m . The 
phase diagram corresponding to this non-competing case is depicted in Fig. 
2(a). 

In the remainder of this paper, we will perform detailed numerical studies of 
two particular aspects of the phase diagrams shown in Figs. 2(a) and 2(b). 
The first concerns the exact location of the multicritical point K m for hi — 0, 
for a finite hs- The second is the verification of the Ising-like nature of the 
low temperature transition lines in Fig. 2(b) for the competing case. 

3 Location of the Multicritical Point 

To quantitatively locate the multicritical point K m , we have performed ex- 
tensive Monte Carlo simulations of Eq. (1), by adapting a modified Wolff 
algorithm. The Wolff algorithm was originally developed for isotropic, con- 
tinuous spin systems such as the XY model ||. In our case, we have added 
symmetry breaking fields to the model. This is accounted for by modify- 
ing the Wolff cluster update algorithm in the following way. We divide the 
Hamiltonian of Eq. ([!]) into two parts, the isotropic part Hxy, 

H XY = -K J2 cosfa-h), (10) 

<i,j> 
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and the anisotropic part i?4,8, 



#4,8 = -/i 4 I>os(40;) + /i 8 ^cos(80,). (11) 

i i 

We form the Wolff cluster for the isotropic part Hxy in the usual fashion us- 
ing a cluster labeling technique similar to the "ants in the labyrinth" scheme 
0. We then calculate the change in the energy due to the anisotropy fields 
for the old and the proposed new cluster as 



Atf 4i8 = - Hfl (12) 

Whether or not the cluster is flipped is determined by applying the standard 
Metropolis acceptance criterion to this energy difference. It is easy to verify 
that this combined algorithm satisfies detailed balance and is ergodic. 



It has recently been shown [|10| that the Wolff algorithm performs poorly 
for anisotropic XY models at low temperatures. It probes the phase space 
effectively but in a presence of strong anisotropy fields, reaching thermal 
equilibrium from an initial state can take very long. Metropolis algorithm, 
on the other hand, reaches local equilibrium rapidly but fails to search the 
phase space extensively. We have overcome these problems by a scheme 
where Wolff and Metropolis algorithms are simply combined by inserting 
several Metropolis local update sweeps after a certain amount of Wolff steps. 



Similar approach has also been suggested previously in Ref. [11 



To locate the multicritical point of the XY model with anisotropy fields, we 



use the method of Lee and Kosterlitz [12]. In this scheme, the "free energy" F 



for a given system of linear size L with order parameter \1/ and with periodic 



10 



boundary conditions can be expressed as 

exp[-F{%L,N)] = NZ- 1 (P)J2^(E,^)^M-PE), (13) 

E 

where N is the number of samples (configurations), Z(/3) is the partition 
function, Q(E, \]/) is the number of states with energy E, and we assume 
that the transition is driven by an external field. In the following, we shall 
for simplicity drop the N dependence of F. F differs from the actual bulk 
free energy but its shape is identical to that of the bulk free energy and 
thus is also the difference AF (see Eq. ([15]) below). The "free energy" has 
two minima due to the two coexisting phases. These minima, located at 
and ^2, are separated by a maximum at \l/ m . In the thermodynamic 
limit the double minima structure vanishes at the transition points but in a 
finite system it may persist even above the transition temperature. In this 
method it is precisely this property that is exploited to reveal the order of 
the transition for a finite system. More specifically, the free energy can be 
expanded below the transition where the correlation length ^ < I as 

F{% L) = L d f (^, g) + L d - l h{m, g) + ..., (14) 

where /o(^, g) is the bulk free energy density, fi(^f, g) is a surface term which 
has a maximum at ^ < \I/ m < \l/ 2 , and g is a scaling field g oc (T — Tc). It 
can then be shown that the free energy has a minimum on both sides of the 
maximum ty m , and that the height difference between the minima and the 
maximum is 
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AF(L) = F(* m , L) - L) = A(g)L d ^ 1 + B{g)L d ~ 2 + ... . (15) 

This expansion holds for first order transitions when £ -C L, and the free 
energy difference AF(L) is an increasing function of the system size L. Even 
for L <C £, AF(L) is an increasing function of L. Thus, at the first order 
transition point when F(^>i, L) = F(ty 2 , L), AF(L) is an increasing function 
of L. In the disordered phase, the free energy difference AF(L) decreases as 
a function of L. 

Near a fixed point describing a continuous transition, a scaling form can be 



developed for the singular part of the free energy [12| . Its analytic expansion 
gives AF(L) = a — bgL 1 ^ + 0(g 2 L 2 ^ u ), where a and b are L-independent 
constants. This form is appropriate for L <C £> where g > for K < K m and 
g < for > i^ m . Thus, we expect that AF(L) increases with L in the 
low temperature phase. However, for our model the behavior in the vicinity 
of K m is expected to be very complicated, and a finite size scaling form has 
not been developed. Naively, we would expect that because 1/u = and 
ln(£) ~ g~ 1//2 , the barrier AF(L) will increase as g\n 2 (L). For K > K m , 
however, AF(L) must increase more rapidly with L eventually crossing over 
to the linear behavior of Eq. (15) for L>( deep into the first order regime. 
We will use this property of AF(L) and change K at /14 = to locate K m , 
as explained below. 

To facilitate the use the finite size scaling technique of Lee and Kosterlitz, 
we calculated the histogram of the order parameter $ = (1/L 2 ) J2i cos40j 
summed over the lattice by dividing the interval [0, 1] into 200 equal bins, 
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and putting each value into its respective bin. Thus, we can construct a 
histogram which shows the double peak structure. This histogram is essen- 
tially an approximation of the partition function. By taking the negative 
of the logarithm of this histogram, we obtain an approximation for the free 
energy distribution F{9,L) of the system (cf. Eq. (|T3|)). As we reside on 
the transition line at = 0, the two peaks of F(^, L) are equally high, and 
we can readily calculate the difference AF(L) given in Eq. (U). For the 
eight fold field, we used the value h$ = —0.15. When using the combined 
algorithm, we first did 3000 Metropolis steps to reach a local equilibrium 
and then continued with 5000 Wolff cluster formations with 10 Metropolis 
steps after each 1000 cluster formations. All this information was discarded. 
The data were averaged over 1 000 000 cluster formations so that after every 
1000 Wolff cluster formations ten Metropolis steps followed. We calculated 
the order parameter histogram for several systems of sizes L = 8, 12, 16, 24, 
32, 48, and 64 at various temperatures along the line /14 = 0. From these data 
we can then deduce AF(L). In the first order regime and in two dimensions, 
it should scale linearly with increasing L. If the transition is continuous, the 
double peak structure should vanish in the limit L — > 00. We also studied the 
distribution of angles by constructing a histogram of each individual angle 

4>i- 

Our main results are depicted in Figs. 3(a)-(c). Typical histograms for 
various systems sizes are shown in Figs. 3(a) and (b), and the extracted 
energy barriers AF(L) as a function of the linear system size L in Fig. 3(c), 
for the present case of non-competing anisotropy fields. All the AF(L)'s were 
calculated by fitting an eight order polynomial to the data. All data points 
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are averages of about 10 configurations. The energy difference increases with 
L for temperatures corresponding to K > 2.3, which indicates a first order 
regime. The behavior of AF(L) at temperatures corresponding to K < 2.2 
indicates, on the other hand, a regime where the transition is continuous. 
At K = 2.2, AF(L) first seems to increase with L up to L = 12 and then 
decrease for larger L, although within error bars it's almost constant. At K — 
2.1, no double peak structure exist. We also analysed the size dependence 
of the multicritical point K m . By fitting AF(L) vs. the logarithm of inverse 
temperature K for L = 8,12,16,24,48,64, we were able to determine the 
multicritical point K m (L) for each system size. From these we estimated the 
multicritical point by scaling this data against 1/L. We conclude that 



^ = 2.1 ±0.1, (16) 

which is the main result of this section. It is also well within the theoretical 
bounds 2/ Ti < K m < 8/n, as expected. 

Finally, we should note that the accuracy of the result suffers from severe fluc- 
tuations in the vicinity of the multicritical point, in particular for the largest 
systems. This inhibits the use of histogram techniques for extrapolation |L3 



4 Ising Transition in the Case of Competing 
Anisotropy Fields 

The other scenario for our model is the case where the anisotropy fields are 
competing. It was shown analytically that the first order transition at low 
temperatures opens up to two Ising transitions with an intermediate phase 
in between. In this intermediate phase the long range order is dictated by 
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the eight fold field. To see this transition numerically, we chose a finite four 
fold field h 4 = 0.06 which favors the orientations fa = 0, 7r/2, etc. for the 
individual spins. The competing eight fold field is chosen to be equal in mag- 
nitude to the one used in the non-competing case i.e. h$ = 0.15. We try to 
locate the corresponding Ising transition temperature Kj l by scanning the 
inverse temperature K. Kj should be well above our estimate of K m « 2.1 
(cf. Fig. 2(b)). For this calculation, we used the Monte Carlo renormal- 
ization group (MCRG) scheme proposed by Binder |14]]. Consider the XY 



spins di = (cos0j, sin0,) on a two-dimensional lattice which is divided into 
subcells or blocks of (linear) size Lb- Let us first define a block variable 

®Lb = JT J2 1>i, (17) 

B i & L B 

where fa is a measure for the local order in the system. We can then define 
an order parameter for each block size as ^ l b =< &l b > where brackets 
again denote a configurational average. We studied different moments of the 
block variables $l s , and constructed the fourth and sixth order cumulants 



Ul b and Vl b for each block size as in Ref. |fL4}| . The variation of these two 
cumulants as a function of the block size Lb gives a flow diagram analogous 
to that of a renormalization group method. These cumulants approach zero 
above Tc as the block size increases. Below Tc, both cumulants tend to 
nonzero values Ul b —> 2/3, and Vl b — > 8/15 as Lb — > 00. At the critical 
point Tc, the cumulants approach nontrivial fixed point values U* and V*. 
Thus, the behavior of the cumulants is reminiscent of the renormalization 
group flows under subsequent transformations of the length scale. One can 
also estimate the correlation length exponent v from the data in the vicinity 
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of U* by noting that 



141 




U* 



U* 




(18) 



for subsystem blocks of size L' B = Lg/b, and by using the scaling relation 
a = 2 — vd. 

For the block variable we chose 



This order parameter, when the angles are folded between — 7r/4 < <fi < 7r/4, 
is zero in the high temperature phase and finite in the low temperature phase 
when a single domain dominates below and above the transition. We studied 
the fourth and sixth order cumulants Ul b an d Vl b • From the flows of these 
cumulants as a function of the inverse linear size of the block, we can deduce 
the nontrivial fixed point value from which we can further extract Kj. 

At low temperatures, the simulations suffer from high barriers between dif- 
ferent regions of the phase space, which we tried to overcome by using the 
combined algorithm. We first used 5000 Metropolis steps to bring the con- 
figuration to a local equilibrium, and then continued with 3000 Wolff cluster 
formations before we started to collect data. The data were averaged over 100 
000 cluster formations, and after every 1000 cluster moves ten thermalizing 
Metropolis steps were completed. 

The flow of the fourth order cumulants is depicted in Fig. 4. We can readily 
see that the value of K at which the transition takes place is Kj e [4.9, 5.0]. 
A more detailed extraction of Ki was done by studying the ratio Uu /Ul b 



1 



B i 



( sin &) • 



(19) 
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where L' B E [L B ,64] and L B = 2,4,6,8,12,16,20,24,32,48,64. The point 
where Ul' b /Ul b = 1 was taken as the transition point. This extraction finally 
gave the result 



Ki = 4.95 ± 0.05. (20) 

By flipping the sign of the four fold field, we also confirmed numerically 
that the Ising transition boundaries are symmetrical with respect to the line 
h 4 = (cf. Fig. 2(b)). It is also possible to use a different order parameter 
(< cos0jsin0j >) by rotating the spins into the first quadrant. The results 
using this order parameter agreed with the choice of < sin0j >. We also 
used the order parameter blocks to estimate the critical exponent v for the 
transition point, and obtained v — 1.1 ± 0.1 in very satisfactory agreement 
with the exact Ising result of v = 1. 

5 Monte Carlo Transfer Matrix Method 

Besides Monte Carlo simulations, another numerical approach which has been 
quite successful in the study of statistical mechanical models is the transfer 



matrix method p!o| . In this method, the free energy of the model defined 
on an infinite strip can be obtained directly from the largest eigenvalue Ao 
of the transfer matrix as — logAo- These calculations can be done exactly 
for models with discrete degrees of freedom when the transfer matrix is of 
low order. When this technique is combined with finite size scaling, one can 
obtain very accurate estimates of critical exponents and other quantities. For 
models with higher order transfer matrix or continuous degrees of freedom, 
as for the case of the XY model with symmetry breaking fields, one has to 
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resort to a Monte Carlo Transfer Matrix (MCTM) technique [16[] in order 



to estimate the largest eigenvalue of the transfer matrix. Some models with 
continuous symmetry have already been studied by this method, including 
frustrated XY [[17], |T^] and coupled XF-Ising models fll9[ . 



We proceed to summarize the method. Further details can found in Refs. [|15 
|T6| , [378] , p0[] . The MCTM consists in a stochastic implementation of the well- 
known power method to obtain the dominant eigenvalue of a matrix. First, 
helical boundary conditions are implemeted in order to get a sparse transfer 
matrix. At each step a new configuration is obtained from the previous 
one, by adding a new spin s' L+l and relabeling the sites. The infinite strip 
can be constructed by repetition of this identical elementary steps. This 
process defines the transfer matrix to add a single site T(s',s), where s = 
s±, S2, sl represent a configuration of a column with L spins. Then a 
sequency of random walkers Ri = si t i, S2,i, ■ ■•Sz.i, 1 < i < r, representing 
the configurations of a column is introduced with corresponding weights W{. 
The number of walkers r is maintained within a few percent of a target 
value r by adjusting the weights properly. A matrix multiplication can be 
regarded as a transition process from s to s' with a probability density P(s', s) 
defined from the elements of the transfer matrix as T(s',s) = D(s)P(s' , s) 
with a normalization factor D(s) independent of s'. In each step the weights 
are changed according to w[ = D(s)wi/c where c = Xor/r is chosen to 
maintain r close to ro, with Ao a running estimate of the eigenvalue. In this 
procedure an MC step consists of a complete sweep over all random walkers. 
After disregarding to MC steps for equilibration, an estimate of the largest 
eigenvalue can be obtained as 



T 

A = , (21) 

t=*o+l 

where = ^ with u^j denoting the configuration weights of a column 
at time t, and T is the total number of MC steps. 

For the calculations of the free energy of our model we performed extensive 
runs using typically r = 20 000 random walkers and T = 100 000 MC steps 
which corresponds to 2 x 10 9 attempts per spin. We concentrated our atten- 
tion in two quantities, the interfacial free energy and the central charge. The 
former can obtained from the free energy per site / = — In A , as calculated 
from the transfer matrix on an infinite strip of width L, by a suitable choice 
of the boundary conditions. If antiperiodic boundary conditions are used 
instead of periodic ones, an interface is favored along the infinite strip and 
the associated interfacial free energy AF(L) can be obtained from the free 
energy difference between periodic and antiperiodic boundary conditions as 

AF(L) = L 2 Af. (22) 

The finite size scaling behavior of this quantity should drastically change 
as a function of temperature depending on the relevance of the symmetry 
breaking fields. At low enough temperatures this interface corresponds to a 
sharp boundary between two different phases due to the presence of the sym- 
metry breaking fields which are relevant at these temperatures, and should 
therefore increase as L d_1 . On the other hand, at higher temperatures when 
the symmetry breaking field is irrelevant, the interface correspond to gradual 
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change, in form of a twist of the angle 0, and AF(L) scales as L d ~ 2 , which 
is roughly a constant at d — 2. In fact, in this regime AF const.) can be 
related to the helicity modulus (divided by k B T) 7 as 



7 = 2AF/vr 2 , (23) 

for large enough L. The change in behavior of the interfacial free energy 
between these regimes can be used to find the multicritical point. 

Another important quantity which can be inferred from the MCTM cal- 
culations, is the central charge c, which classifies the possible conformally 
invariant critical theories [ZIJ. For example, c = 1/2 for the Ising model and 



c = 1 along the critical line of the XY model. This quantity can be related 
to the amplitude of the singular part of the free energy per site (at criticality) 
of the infinite strip by 



f(K c ,L)=f + — (24) 

for sufficiently large L, where fo denotes the regular contribution to the free 
energy. Although c is only defined at criticality, the value of c extracted from 
a f(K, L) x 1/L 2 fitting of the free energy as a function of system size, can be 
used to define an effective size and coupling dependent central charge c(K, L) 
away from the critical point. According to the Zamolodchikov c-theorem [[22|j . 
this quantity should have a well defined behavior near the critical point and 
reach a constant value, equal to the central charge, at the fixed point. As 
a consequence c(K, L) should have a maximum near an unstable fixed point 
and converge to c = in the completely disordered or ordered phases. This 
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property is particularly useful in locating the multicritical point and the Ising 
transition. 

To facilitate direct comparison with the MC results of Sec. 3, we set h± = 
and h 8 = —0.15. The results of the MCTM calculations for the central charge 
in this non-competing case (cf. Fig. 2(a)) are shown in Fig. 5. The value 
of c was estimated by fitting the free energy per site to Eq. (24), using sizes 

L = 5 11. First, the abrupt onset of c near K « 1 agrees well with the 

known result for the XY transition (where h$ is irrelevant). The multicritical 
point, on the other hand is estimated to be at the point where the central 
charge begins to decrease from the value of one (the XY phase value) to zero 
at K m = 2.6 ± 0.4. This is in fair agreement with the result of Eq. (16) as 
obtained from the histogram method; the discrepancy can be attributed to 
severe finite size effects for the relatively small strip widths studied. 

In Fig. 6, we show the results for the interfacial free energy for different strip 
widths L as a function of temperature for the same non-competing case. One 
clearly sees a size dependence at large values of K indicating the relevance 
of the symmetry breaking field h 8 . At temperatures above K~} 1/2.6 the 
size dependence is practically absent indicating the spin wave regime of the 
critical line of the XY model. In this regime the helicity modulus is equal 
to the renormalized coupling constant K* and should be equal to K* = 8/n 
at K m . As indicated in Fig. 6, the helicity modulus 7 = 2AF/ir 2 is in good 
agreement with the expected result. Moreover, the termination point of the 
lines at AF = is in good agreement with the KT transition point. 

Next, we use the MCTM method to study the Ising transition in the case of 
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competing h± and h$ fields. As in Sec. 4, we set h± = 0.06 and hg = 0.15. 
The results for c vs. K are shown in Fig. 7. The effective central charge 
has two peaks as a function of temperature. At high temperatures the peak 
value is very close to c = 1, which is consistent with a transition in the 
universality class of the XY model with a four fold field. The other peak 
at low temperatures, is close to c = 1/2 which is consistent with an Ising 
transition. The Ising transition can be estimated simply from the peak value 
of c yielding Kj = 4.6 ±0.4. These results are again in reasonable agreement 
with the MCRG method, Eq. (20). 

6 Summary and Discussion 

In this work, we have analysed in detail the properties of the two dimensional 
XY model in the presence of four fold and eight fold symmetry breaking fields 
hi and h$, respectively. First we have applied mean field and renormalization 
group arguments to predict that when the field changes from positive to 
negative values, there is a phase transition between a phase with the order 
parameter pointing along the <fi = 7r/4 direction, and another phase with 
the corresponding order parameter pointing along the = direction. This 
phase transition is continuous at high temperatures where the eight fold 
field is irrelevant. In this case, right at the phase boundary all the symmetry 
breaking fields are absent and the system is in the ordered phase of a pure XY 
model with only algebraic long range order. At lower temperatures when the 
eight fold field is relevant, the nature of the transition now depends crucially 
on whether h$ and h± are compatible or competing with each other. In the 
former case, the transition becomes first order and there is a multicritical 
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temperature separating the high temperature and low temperature phase 
boundary. In the competing case the first order phase transition splits into 
two transitions both of which are in the Ising universality class. The order 
parameter rotates continuously from the = direction towards the = 7r/4 
direction in between the two Ising phase boundaries. We note that similar 
effects can result from the competition of six fold and twelve fold fields as 
discussed by Selinger and Nelson J5J in their study of liquid crystals. 

We have then performed detailed numerical analysis of our model Hamilto- 
nian through Monte Carlo simulations. Because of the large finite size effects, 
the finite size simulation data has to be extrapolated to infinite size through 
finite size scaling concepts and numerical renormalization group analysis. We 
have first confirmed the qualitative nature of the phase boundary as discussed 
above. We have then employed the recently developed histogram techniques 
to locate the multicritical temperature. Finally, in the case of the compet- 
ing fields, we have used Binder's block renomalization technique to identify 
the Ising like transitions. These results have further been corroborated by 
Monte Carlo Transfer Matrix techniques. All our numerical results are in 
good agreement with theoretical arguments. 

Our model was originally motivated by the study of the H/W(100) chemisorp- 
tion system. We have shown previously that the critical properties of this 
system in the ordered c(2 x 2) phase can be desribed by an XY model 
with four and eight fold symmetry breaking fields. The effect of the ad- 
sorbed hydrogen in the low coverage limit is to change the effective four fold 
field. Thus increasing the hydrogen coverage in this system is tantamount 
to changing the four fold field from negative to positive values in the model 
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Hamiltonian studied in this paper. In addition, our previous work showed 
that the adsorbed hydrogen also generates an eight fold field which is 
compatible with the four fold field. Experimental evidence from an infrared 
spectroscopy study for this system supports the scenario of a continuous 
transition from < 11 > (corresponding to negative h±) to < 10 > (positive 
phase at room temperature when the hydrogen coverage is increased. On 
the other hand, the Low Energy Electron Diffraction study of Griffiths et al. 
24fl showed indications of coexistence between the two phases for coverages 



in the range from about 0.05 to 0.16 which indicates a first order transition. 
This is exactly the behavior of our model Hamiltonian studied in this paper 
when the fields /14 and h$ are compatible. Thus the qualitative agreement 
between the experimental data and the theory presented here is gratifying. 
In view of the very rich behaviour of the phase diagram, especially near the 
multicritical point, more experimental studies of the switching transition in 
this system and comparison with the theoretical predictions here would be 
most fruitful. 
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Figure Captions 



Fig. l.(a) A schematic figure of the anisotropy potential V(4>) of Eq. (2), for 
the case of non-competing fields with h 8 = — 1. The curves from bottom to 
the top are for = —6, —4, —2, 0, 2, 4, and 6, respectively. The transition 
at hi = is abrupt, (b) The anisotropy potential for competing fields with 
h$ = +1. The curves are for the same values of as in (a). The transition 
around — is continuous. 

Fig. 2. (a) A schematic phase diagram for the Hamiltonian Eq. (1) on the 
h 4 — K~ x plane, with a fixed non-competing h 8 < 0. For finite values of 
h^, the two transition lines beyond the Kosterlitz-Thouless transition point 
Krt < 7r /2 belong to the universality class of the XY model with a four 
fold symmetry breaking field. Between K^ T and the multicritical point K^ 1 , 
there is a line of pure XY transitions. Below K^ 1 , the eight fold field is 
relevant and there is a line of first order transitions. The lower bound for 
K^ 1 is 7r/8 (cf. Eq. (9)). (b) The counterpart of the phase diagram of (a) for 
the case of a fixed, competing h 8 > 0. Below -ft'" 1 , where h 8 is relevant, the 
first order transitions open up into two lines of continuous Ising transitions, 
which terminate at ±4|/i 8 |. 

Fig. 3. (a) Logarithm of the histogram of the local order parameter $ = 
(l/^ 2 ) Ei cos40i (cf. Eq. (|1D), at K = 2.1 for system sizes L = 8, 12, 16, 48, 64, 
with /i 4 = 0.06 and h s = —0.15. iV is the number of samples in each bin. (b) 
Logarithm of the histogram of the local order parameter $ = (1/L) X)« cos40j 
at K = 2.4 for system sizes L = 8,12,16,24,48,64, with h 4 = 0.06 and 
h s = — 0.15. (c) The free energy barriers AF(L) (cf. Eq. (15)) extracted from 
polynomial fits to the histograms of the order parameter, plotted against L. 
Here = and h 8 = —0.15, corresponding to Fig. 2(a). Above K m w 2.2, 
the barriers clearly start increasing indicating the onset of a first order tran- 
sition. See text for details. 
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Fig. 4. Typical MCRG flows of the fourth order cumulant vs. 1/L for the 
competing case, with h± = 0.06 and h 8 = 0.15. When K is varied, the Ising 
transition occurs at around K I 4.9 — 5.0 (cf. Fig. 2(b)). See text for 
details. 

Fig. 5. MCTM results for the effective central charge along the line h± = 0, 
with h 8 = —0.15 (the non-competing case of Fig. 2(a)). The first abrupt 
onset to c = 1 corresponds to the KT transition, and the onset of decline of 
c at K m ~ 2.6 indicates entry to the line of first order phase transitions just 
below K' 1 in Fig. 2(a). 

Fig. 6. The interfacial free energy from the MCTM method for different 
strip sizes as a function of K, for the case h 4 = and h 8 = —0.15 as in 
Fig. 5. Below K m ps 2.6, the lack of size dependence of the energy indicates 
the onset the spin wave regime corresponding to a line of XY transitions. 
The corresponding value of K* is in good agreement with the theoretical 
prediction 8/n. See text for details. 

Fig. 7. MCTM results for the effective central charge with h A = 0.06, 
h 8 = 0.15, corresponding to the competing case of Fig. 2(b). The first peak 
at c = 1 corresponds to the KT transition point, while the second peak at 
Kj p^ 4.6 verifies the expected Ising transition, with c = 1/2. 
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